## Replication data for "Measuring the Consolidation of Power in Non-Democracies"
## Jennifer Gandhi and Jane L. Sumner
## Journal of Politics
## Date: 4/29/2019
## Please contact Jane Sumner (jlsumner@umn.edu) with questions about replication data and code.


## Load in required libraries
library(reshape2)
library(R2jags)
library(stringr)
setwd("D:/Dropbox/Personalism/Data")

## Load in data
mod.dat <- read.csv("gandhi-sumner-replication.csv",stringsAsFactors=F)

## Model
# Reformat data
st <- melt(mod.dat,id.vars=names(mod.dat[,c(1:5)]),
           measure.vars=names(mod.dat[,6:ncol(mod.dat)]))


# Just in case, making sure they're in the right temporal order.
st.new <- NULL
for(i in 1:length(unique(st$country))){
  st.small <- st[which(st$country==unique(st$country)[i]),]
  st.small2 <- st.small[order(st.small$year),]
  st.new <- rbind(st.new,st.small2)
}
st <- st.new
st <- st[!is.na(st$value),] # remove NAs


N <- nrow(st)  # number of observations in the model
J <- length(unique(st$variable)) # number of observable indicators
I <- length(unique(st$countryleaderyear)) # number of country-leader-years
itemid <- as.numeric(as.factor(st$variable)) # unique indicator id
y <- st$value # observable behaviors
clyid <- NULL
for(i in 1:I){ # indexes unique country-year-regime
  howmany <- nrow(st[which(st$countryleaderyear==unique(st$countryleaderyear)[i]),])
  clyid <- c(clyid,rep(i,howmany))
}
K <- max(clyid)
st$cont <- ifelse(st$firstyear==1,0,1) # indicates first year of regime (for priors)
cont <- unique(st[,c("countryleaderyear","cont")])$cont # indicates continuation of regime (not first year)

irt.ord.model <- function(){
  for (i in 1:N) { # observations
    
    
    y[i] ~ dbern(eta[i])
    logit(eta[i]) <- alpha[itemid[i]] * (theta[clyid[i]] - beta[itemid[i]])
    # }
    
    
  }
  theta[1] ~ dnorm(0,1)  # first observation
  
  sdcont <- 10 # tighter prior for continuing regime
  sdfy <- 1    # looser prior for first-year of regime
  
    for(k in 2:K){
    # is this a continuing regime? if yes == 1, if no == 0
    # if 1, dnorm(x[t-1],.1)
    # if 0, dorm(0,1)
    theta[k] ~ dnorm(theta[k-1] * cont[k],cont[k]*sdcont + (1-cont[k])*sdfy);T(-3,3)
  }
  

  
  for (j in 1:J) { # items
    beta[j] ~ dnorm(0,10) # difficulty
    alpha[j] ~ dnorm(0,10);T(0,) # discrimination
  }     
  
}


inits <-  function(){ list("theta"=runif(K,-3,3),
                           "beta"=runif(J,-1,1),
                           "alpha"=runif(J,.05,1)
)}

jags_data <- list(y=y,clyid=clyid,itemid=itemid,cont=cont,N=N,J=J,K=K)


total.iter <- 50000
jags.mod <- jags(model.file = irt.ord.model,
                     data=names(jags_data),
                     n.chains=5,
                     inits=inits,
                     n.iter=total.iter,
                     n.burnin=total.iter*.95,
                     parameters.to.save=c("beta","alpha","theta"))

end.time <- Sys.time()

# checking Rhat
sum(jags.mod$BUGSoutput$summary[,"Rhat"]<=1.01)/length(jags.mod$BUGSoutput$summary[,"Rhat"])*100
hist(jags.mod$BUGSoutput$summary[,"Rhat"])

# Extracting results and merging with data
out <- jags.mod.C23$BUGSoutput

results <- data.frame(id=unique(clyid),
                      countryleaderyear=unique(as.character(st$countryleaderyear)))
results$year <- unlist(lapply(results$countryleaderyear, function(x){
  unlist(str_split(x,
                   "-"))[length(unlist(str_split(x,"-")))]
}))

results$xhatmean <- apply(out$sims.list$theta,2,mean)
results$xhatmedian <- apply(out$sims.list$theta,2,median)
results$xhatsd <- apply(out$sims.list$theta,2,sd)
results$xhat025 <- apply(out$sims.list$theta,2,function(z){quantile(z,.025)})
results$xhat975 <- apply(out$sims.list$theta,2,function(z){quantile(z,.975)})
# 

mod.dat$countryleaderyear <- paste(mod.dat$leaderspell,mod.dat$year,sep="-")
mod.dat2 <- merge(mod.dat,results,by=c("countryleaderyear","year"))
results <- mod.dat2

## Figure 1: Difficulty and Discrimination Plots

# if (and only if) you have estimated the model above, run this code as well:
item.chars <- data.frame(item=names(mod.dat)[6:ncol(mod.dat)],
                         discrimination=apply(out$sims.list$alpha,2,mean),
                         discrim025=apply(out$sims.list$alpha,2,function(z){quantile(z,.025)}),
                         discrim975=apply(out$sims.list$alpha,2,function(z){quantile(z,.975)}),
                         difficulty=apply(out$sims.list$beta,2,mean),
                         diff025=apply(out$sims.list$beta,2,function(z){quantile(z,.025)}),
                         diff975=apply(out$sims.list$beta,2,function(z){quantile(z,.975)}))


# if you have not estimated the model and want to run the plots with our results,
# run these lines of code to load in our results from the paper:
results <- read.csv("results-1-4-2019.csv",stringsAsFactors=F)
item.chars <- read.csv("item-chars-1-4-2019.csv",stringsAsFactors=F)


item.chars$name <- c("Domestic purges","Military purges","Cabinet change",
                     "Dictator before transition",
                     "No formal collective",
                     "Military involvement",
                     "No regime party",
                     "Founded while in power",
                     "Multiple parties",
                     "No position",
                     "Two positions",
                     "3+ positions",
                     "Family in power")

#png("difficulty.png")
par(mar=c(4,10,2,1),mgp=c(2,1,0))
plot(item.chars$difficulty,c(1:nrow(item.chars)),pch=16,
     yaxt="n",ylab="",xlab="estimate",main="Difficulty Parameters",
     xlim=c(min(item.chars$diff025),max(item.chars$diff975)))
segments(y0=c(1:nrow(item.chars)),x0=item.chars$diff025,x1=item.chars$diff975)
axis(2,at=c(1:nrow(item.chars)),item.chars$name,las=2,cex.axis=.8)
abline(v=0,lty=3)
#dev.off()

#png("discrimination.png")
par(mar=c(4,10,2,1))
plot(item.chars$discrimination,c(1:nrow(item.chars)),pch=16,
     yaxt="n",ylab="",xlab="estimate",main="Discrimination Parameters",
     xlim=c(min(item.chars$discrim025),max(item.chars$discrim975)))
segments(y0=c(1:nrow(item.chars)),x0=item.chars$discrim025,x1=item.chars$discrim975)
axis(2,at=c(1:nrow(item.chars)),item.chars$name,las=2,cex.axis=.8)
abline(v=0,lty=3)
#dev.off()

## Figure 2: Plots for Ethiopia and Tanzania
# This function will make a plot for any country in the data.
mod.dat$leadername <- str_trim(gsub("[(].*","",mod.dat$leaderspell))
mod.dat$leadercountry <- paste(mod.dat$leadername," (",mod.dat$country,")",sep="")
countryplot <- function(countryname,name=countryname){
  par(mar=c(4.1,2.1,2.1,3.1))
  res.small <- results[which(results$country==countryname),]
  cols <- "black"
  plot(res.small$year,res.small$xhatmean,pch=16,
       main=name,ylab="personalism estimate",xlab="year",col=cols,
       ylim=c(min(results$xhat025)-.5,max(results$xhat975)+.5))
  abline(v=(as.numeric(res.small$year[which(res.small$firstyear==1)])-.5),
         lty=3)
  res.small$leadername <- gsub(" [(].*","",res.small$leadercountry)
  for(j in 1:length(unique(res.small$leadername))){
    yrs <- as.numeric(res.small$year[which(res.small$leadername==
                                               unique(res.small$leadername)[j])])
    if(length(yrs)>3){
      text(min(yrs)+(max(yrs)-min(yrs))/2,
           (max(res.small$xhat975[which(res.small$leadername==
                                          unique(res.small$leadername)[j])])+.1),
           gsub(" ","\n",unique(res.small$leadername)[j]),cex=.8)
    }
  }
  
  segments(x0=as.numeric(res.small$year),
           y0=res.small$xhat025,y1=res.small$xhat975,col=cols)
  
}

par(mfrow=c(1,2))
countryplot("Ethiopia")
countryplot("Tanzania, United Republic of","Tanzania")
par(mfrow=c(1,1))

## Figure 3: Cross-section variance in power consolidation estimates during the 1946-2008 period
# Variation cross-sectional v temporal?
plot(c(1946:2008),tapply(results$xhatmean,results$year,var),pch=16,type="h",lwd=2,
     main="",ylab="Variance",xlab="Year")
tapply(results$xhatmean,results$year.x,var)[order(tapply(results$xhatmean,results$year.x,var))]
cor(tapply(results$xhatmean,results$year.x,var),tapply(results$xhatmean,results$year.x,length))


## All tables calculated in STATA. See associated STATA file.

## end of file.